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Chapter  1 


Overview 


1.1  Identification  of  the  Problem 

After  several  decades  of  work,  many  unresolved  issues  remain  associated  with  the  generation 
and  control  of  thermoacoustic  combustion  instabilities  in  full-scale  combustors.  Academia 
has  been  working  on  theoretical,  numerical,  and  experimental  aspects  of  combustion  insta- 
biltiy  which  can  be  transferred  to  immediate  relevance  for  design  engineers  in  the  combustion 
industry.  From  the  industrial  perspective,  an  urgent  need  exists  for  the  development  of  us¬ 
able  models  which  can  predict,  at  some  defined  level  of  certainty,  the  potential  occurrence 
of  combustion  instabilities  over  the  entire  operating  range  of  input  conditions.  Because  of 
the  nonlinear,  coupled,  and  very  complex  nature  of  the  physical  processes  underlying  the 
instabilities,  it  is  highly  desirable  that  these  models  be  developed  using  simplified,  reduced- 
order  structures.  Modeling  efforts  of  this  type  are  ongoing.  However,  the  complexity  of  the 
unstable  thermoacoustic  dynamics  may  also  be  addressed  through  the  use  of  more  detailed 
CFD  modeling  approaches.  This  approach  has  the  distinct  advantage  of  providing  field 
information  which  is  especially  important  in  active  combustion  control  (ACC). 

Advances  in  CFD  software  and  related  hardware  have  greatly  enhanced  the  prospects 
for  physies-bascd  modeling  of,  cb nm: sally  ^reacting,  flows  .*• Naturally,  .as  the  CFD ,  field  su¬ 
tures,  its  influence  on  design  and  analysis  in  advanced  aerospace  applications  becomes  more 
prevalent.  For  example,  the  structured  flow  package,  QASV^  licensed  by  AeroSoft,  provides 
a  wide  variety  of  user-controlled  options  for  modeling  subsonic,  transonic,  supersonic  and 
hypersonic  flows.  These  include  an  extensive  chemistry  database  containing  species,  reac¬ 
tions  and  models  of  interest,  options  for  frozen  or  non-equilibrium  chemistry,  and  multiple 
options  for  thermodynamics  and  turbulence  modeling. 

From  experience  with  augmentors,  we  know  that  bluff-body  burners  equipped  with  mul¬ 
tiple  fuel  injectors  can  increase  fuel/air  mixing  and,  therefore,  could  be  used  effectively  in 
the  main  burner  as  well.  However,  thermoacoustic  instabilities  can  cause  structural  dam¬ 
age  when  acoustic  pressure  fluctuations  couple  with  unsteady  heat  release  from  combustion. 
Light-weight  active  control  systems,  which  focus  on  the  mixing  behind  these  bluff-body  burn¬ 
ers,  can  alter  the  feedback  mechanism  and  reduce  the  thermo-acoustic  instability.  Defining 
such  a  system  requires  a  tool  for  determining  generic  combustor  characteristics  which  inspires 
for  our  approach. 

We  make  a  concerted  effort  to  develop  and  test  an  unsteady,  sensitivity  tool  using  high- 
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■  CHAPTER  1.  OVERVIEW 


Figure  1.1:  Fuel  nozzle  showing  the  range  of  positions  for  the  swirling  vanes. 


order  computational  simulations,  experimental  results  for  verification/validation,  and  opti¬ 
mal  control  laws  to  control  combustion  flow  fields  susceptible  to  thermo-acoustic  instabilities. 


1.2  Technical  Problem  Addressed  in  Phase  I 


The  Phase  I  effort  addressed  the  computational  problem  of  modeling  steady  and  unsteady 
combustion  in  a  generic  dump  combustor  with  swirl  (see  Figure  1.1).  We  studied  two  com¬ 
ponents  of  the  combustor  computationally:  (1)  the  swirling  air/methane  mixing  inside  the 
fuel  nozzle,  and  (2)  the  response  in  the  combustion  region  to  a  time-dependent  velocity 
profile.  We  were  interested  in  the  sensitivity  of  the  flow  field  with  respect  to  the  axial  posi¬ 
tion  of  the  swirling  vanes  and  to  the  axial  jet  velocity.  We  solved  the  flow  equations  using 
AeroSoft’s  QASV  flow  solver  and  the  sensitivity  equations  using  the  Sensitivity  Equation 
Method  (SEM)  incorporated  into  AeroSoft’s  SSAfSS  software.  We  identified  an  effectivean 
elementary  chemical-kinetics  model  which  uses  32  reactions  and  12  species.  The  project, 
however,  did  not  solve  the  unsteady  sensitivity  equations,  which  inspires  part  of  the  Phase  II 
plan. 


1.3.  RELATIONSHIP  TO  PHASE  II  . . : .  ‘  "  '  "  •" . 3 

1.3  Relationship  to  Phase  II 

The  selection  of  a  chemistry  model  and  the  implementation  of  that  model  into  the  QA.SV 
flow  solver  and  SENSE  sensitivity  code  was  the  main  accomplishment  of  the  Phase  I.  We 
uncovered  several  computational  issues  to  improve  including  tight-coupling  of  the  sensitivity 
solver  into  the  flow  solver  for  true  unsteady  sensitivities  and  the  incorporation  of  methods 
to  relieve  the  stiffness  of  chemically  reacting,  time-dependent  problems.  Additionally,  a 
coupling  mechanism  must  exist  between  acoustic  and  heat-release  perturbations  to  accurately 
model  combustion  instabilities.  Through  experience  and  research  of  related  combustion 
literature,  we  emphasize  the  importance  of  internal-flow  boundary  conditions  upon  the  self¬ 
excitation  of  the  instability.  These  are  issues  which  we  address  in  Phase  II. 


1.4  Project  Objectives 


The  overall  goal  of  this  project  is  to  provide  a  sensitivity-analysis  tool  for  the  control  of 
heat-release  rate  distribution  in  aeroengine  combustors  with  emphasis  on  bluff-body  type 
flame-holders.  A  bluff  body  serves  as  the  anchoring  mechanism  and  we  are  interested  in  the 
best  way  to  introduce  fuel  into  the  combusting  flow.  This  report  documents  the  objectives, 
work  carried  out  and  results  obtained  during  the  Phase  II  project.  To  achieve  the  Phase  II 
goals,  the  following  technical  objectives  are  given  below: 


1.  Using  parametric,  scaled  experiments  on  a  2-D,  bluff-body  combustor,  we  obtain  a 
complete  characterization  of  the  steady  and  unsteady  conditions.  Additionally,  a  para¬ 
metric  study  of  the  combustion  performance  will  be  conducted  under  open-loop  control. 
Data  from  these  experiments  will  help  test  the  computational  tool. 


2.  Implement  a  variety  of  computational  elements  for  simulating  combustion  instabili¬ 
ties:  (1)  time-accurate  numerics  which  address  chemically  reacting  flows,  (2)  a  two- 
parameter  coupling  mechanism  between  acoustic  and  convective  disturbances,  and 
(3)  non-reflective  boundary  conditions  to  remove  numerical  noise. 


3. 


Combine  developed  seifSitivity  technologies'  withPai  cOEftrnerci&i'fioW^sdiVer'  package  to*  - 
compute  sensitivities  to  actuator  parameters.  Truly  unsteady  sensitivities  can  only  be 
obtained  by  tightly  coupling  the  flow  and  sensitivity  solvers.  Therefore,  we  couple  the 
two  and  beneficially  enhance  the  commercial  appeal  of  QASV  in  the  process. 


j.vfc.yu. L' 


4.  A  distributed-parameter  control  model  has  been  formulated  to  qualitatively  capture 
thermo-acoustic  instability  features  in  a  combustion  process.  Tools  for  numerical  sim¬ 
ulation  and  stability  analysis  have  been  developed  and  demonstrated.  Feedback  stabi¬ 
lization  has  been  investigated  through  two  optimal  control  problems. 


CHAPTER  1.  OVERVIEW 


*  1.7.T 


Chapter  2 

Phase  II  Work  and  Reults 


2.1  Experimental  Study  of  Bluff-Body  Combustor 


Conceptual  Design 


In  the  case  of  bluff-body  stabilized  combustion,  strong  coupling  between  acoustic  and  shear- 
layer  instabilities  is  expected.  This  will  manifest  itself  in  the  shedding  of  large-scale  struc¬ 
tures,  which  are  typically  straddled  by  the  flame/combustion  zone.  In  the  conceptual  stage, 
we  used  simple,  cold-flow  CFD  simulations  to  obtain  an  approximation  of  these  flow  struc¬ 
tures.  The  geometrical  configuration  was  that  of  the  coaxial  bluff-body  combustor  shown 
schematically  in  Figure  2.1.  This  calculation  assumes  two-dimensional  flow  and  is  modeled 
using  the  RNG  K-e  model  with  second-order  accuracy  in  time  and  space. 

The  most  important  piece  of  information  from  the  calculation  is  the  frequency  content 
of  the  flow.  A  time  series  of  vorticity  magnitude  is  shown  in  Figure  2.2.  We  can  clearly 
see  the  shedding  of  alternating  vortices.  The  frequency  contents  of  the  turbulent  flow- 
field  was  probed  (from  numerical  data)  at  six  locations  downstream  of  the  bluff  body.  A 
magnification  of  the  spectra  is  shown  in  Figure  2.3  which  indicates  an  existence  of  both 
odd  and  even  harmonics  for  the  lower  inlet- velocity  case  ( i.e 15 m/s).  For  the  higher 
•velocity ’case  only 4he;evemliariiionks 30  =• 

the  fundamental  vortex  shedding  frequency  corresponds  to  a  Strouhal  number  of  0.3.  Note 
that  from  Figure  2.2,  the  time  taken  for  one  vortex  to  shed  is  approximately  8  x  10“3  s, 
which  corresponds  to  a  frequency  of  125  Hz.  The  frequency  calculated  by  FFT  is  120  Hz 
(for  the  15  m/s  case)  which  corresponds  to  a  shedding  time  of  approximately  8.33  x  10“3  s. 


Preliminary  Design 

With  the  cold-flow  calculations  complete,  we  moved  to  design  a  combustor  for  acoustic 
frequencies  in  the  range  of  50-500  Hz.  Because  we  wish  to  excite  different  frequencies, 
the  design  calls  for  an  ability  to  change  the  combustor’s  acoustic  characteristics  by  either 
changing  the  bluff-body  mounting  arm  or  adding  downstream  sections  to  the  test  facility. 
Initial  concepts  of  mounting  the  bluff-body  on  an  arm  guided  by  slip  bearings  was  discounted 
because  of  cost  and  construction  complications.  Maximum  optical  access  of  the  test  section 
is  desired  which  includes  ports  for  upstream  and  downstream  acoustic  excitation. 

The  preliminary  design  of  the  test  rig  is  shown  in  Figure  2.4.  The  rig  is  fully  modular 
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Figure  2.1:  Coaxial  bluff-body  combustor  geometry  used  in  the  CFD  simulation.  Dimensions: 
D  =  76.2  mm.,  Lp  =  3D,  d  =  D/2. 


and  comprises  flow  conditioning,  a  combustor  section  and  an  exhaust  mechanism.  The 
ng  can  be  operated  in  both  pre- 1  nixed  and  di ffu  sions-Tho de? .  The  combust : on  .section  har 
a  circular  geometry  with  a  diameter  of  101.6  mm  (4  in)  and  a  length  of  508  mm.  The 
flame  section  is  air-cooled  and  double- walled.  Quartz  windows  are  provided  on  three  sides 
for  viewing-access  and  optical  diagnosis.  The  chamber  can  be  acoustically  modified  by  the 
addition  of  downstream  sections.  A  back-pressure  valve  in  the  exhaust  section  provides  the 
capability  of  operating  the  test  rig  at  elevated  pressures.  The  water-cooled  bluff  body  is 
shown  in  Figure  2.5.  The  disc  measures  50.8  mm  in  diameter  and  the  body  if  fitted  with 
one  axial-injection  hole  and  four  radial-injection  holes,  each  located  90°  apart.  Each  of  the 
five  jets  has  an  independent  fuel  supply  and  control. 

A  second  CFD  investigation  was  performed  for  the  preliminary  combustor  under  hot-flow 
conditions.  Parametric  studies  were  conducted  with  mean  combustor  inlet  velocities  of  30  and 
15  m/s  to  determine  the  size  and  shape  of  the  bluff  body  and  size  of  the  injection  holes.  From 
the  initial  cold-flow  study,  an  incoming  inlet  velocity  of  15  m/s  vortex-shedding  frequency 
of  120  Hz  (corresponding  to  a  Strouhal  number  of  0.3).  Reacting  flow  was  simulated  for 
the  geometry  shown  in  Figure  2.6.  Three  injection  holes  were  included  in  the  bluff  body 
geometry  to  account  for  the  methane  fuel  injection  from  the  bluff  body.  Two  radial  holes  on 
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t  =  8  x  10-3 s  t  =  9  x  10  3  s 


Figure  2.2:  Vorticity  magnitude  contours  ( Uin\et  =  15  m/s) 
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Frequency  (Hz) 


Figure  2.3:  Power  spectrum  plot  of  vorticity  magnitude  (Ptll;  Uiniet  =  15  m/s  and  30  m/s). 


Figure  2.4:  Schematic  of  the  combustor  experimental  set-up. 
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Figure  2.5:  Schematic  of  the  water-cooled  bluff  body. 


the  top/bottom  of  the  bluff  body  and  one  axial  hole  are  considered. 

A  large  eddy  simulation  (LES)  of  the  reacting  flow  was  carried  out  using  the  localized 
dynamic  sub-grid  model  (LDKM)  and  a  15-step  reduced  reaction  mechanism  for  methane-air 
premixed  combustion  (where  the  equivalence  ratio  was  <j)  =  0.5).  The  calculations  predict 
that  the  flame  in  anchored  by  the  bluff  body  and  that  no  vortex  shedding  was  present. 
Even  though  combustion  tends  to  homogenize  the  vortex  shedding  process,  these  results 
were  considered  preliminary  and  therefore  a  final  conclusion  could  not.be  dr^wn  concerning 
vortex  sb“dding.  • 

With  preliminary  design  work  and  reacting-flow  simulations  complete,  efforts  then  fo¬ 
cused  on  the  development  of  an  experimental  combustion  test  facility  at  Virginia  Tech.  A 


17  in 


Figure  2.6:  Internal  geometry  for  the  CFD  analysis 
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Mean  Inlet  Cross-flow  Velocity  of  Premixed  Charge  (m/sj 

15 

Equivalence  Ratio 

0.5 

Cross  Flow  Density  (kg/m3) 

1.155 

Air  Mass  Flow  Rate  (kg/s) 

0.138 

Air  Volume  Flow  Rate  (scfm) 

250.458 

Fuel  Mass  Flow  Rate  (kg/s) 

0.004 

Fuel  Volume  Flow  Rate  (scfm) 

13.204 

Total  Volume  Flow  Rate  in  Pre-mixer  (scfm) 

263.663 

Cross  Flow  Reynolds  No 

66671.481 

Cross  flow  Velocity  at  Radial  Injection  Plane  (m/s) 

19.914 

Injection  Diameter  (mm) 

2 

Modulation  Mass  Flow  Rate  (kg/s) 

10%  of  the  fuel  in  the  cross  flow 

Modulation  Volume  Flow  Rate  (scfm) 

1.320 

Modulation  Jet  Velocity  (m/s) 

126.936 

Modulation  Jet  Reynolds  No 

18865.844 

Thermal  Power  of  the  Rig  (W) 

222117.188 

Table  2.1:  Combustor  rig  parametrics. 


facility  capable  of  providing  1200  scfm  of  vitiated  airflow  at  190  psig  (which  can  be  enhanced 
to  35  bar  maximum)  pressure  and  1050° F  temperature  was  established.  On  the  fuel  side, 
natural  gas  supplies  at  750  psig  was  constructed.  A  desire  to  automate  the  operation  of 
the  facility  in  terms  of  supervision  and  running  the  rig  was  completed.  Detailed  engineering 
drawings  of  the  combustor  test  rig  (portions  of  which  are  shown  in  Figure  2.7(a)  and  Fig¬ 
ure  2.7(b))  were  completed  and  sent  to  Virginia  Tech’s  machining  operation  for  milling  and 
fabrication. 


Design  Rig  Summary  The  rig  is  187  in.  long  and  can  be  operated  at  5  bar  maximum 


provided  i 

acoustically  forced  or  modified  to  be  self-excited.  The  water-cooled  bluff  body  has  a  disc 
diameter  of  2  in.  with  one  axial  and  four  (90°  apart)  radial  injection  holes,  each  with  inde¬ 
pendent  fuel  supply  and  control.  The  main  operating  parameters  of  the  combustor  are  given 
in  Table  2.1. 


Combustor  Fabrication 

The  fabrication  of  the  high-pressure  combustor  rig  took  place  in  two  locations.  The  ma¬ 
jority  of  the  parts  were  fabricated  and  assembled  at  the  Virginia  Tech  Mechanical  Engi¬ 
neering  (ME)  Department’s  Machine  Shop.  The  combustor  balance  was  constructed  by  a 
University  contractor. 

At  this  point,  the  experimental  portion  of  the  project  began  experiencing  set-backs. 
When  the  rig  became  65%  complete,  budget  cuts  by  the  state  of  Virginia  began  to  hamper 
progress.  The  ME  portion  of  the  job  was  done  free  of  charge  to  the  project.  However,  budget 
costs  constrained  the  amount  of  man  power  and  time  that  was  assigned  to  the  project.  Once 
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(a)  3-D  view  of  rig. 


- 187.1975 - 

(b)  Dimensions  of  rig. 


Figure  2.7:  Bluff-body  combustor  view  and  dimensions. 
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the  machining  activities  ended,  welding  of  the  large  stainless  steel  flanges  commenced  by  an 
ASME-certified  outfit  in  Salem,  Virginia.  Some  of  the  parts  are  shown  in  Figure  2.8  and 
Figure  2.9. 

After  many  delays,  assembly  of  the  combustor  rig  was  accomplished  (see  Figure  2.10  for 
the  assembled  combustor  rig).  The  compressed  air,  fuel  and  water  (for  cooling  the  bluff  body 
and  the  exhaust  gases)  feed  systems  were  checked  for  high-pressure  leaks  and  functionality 
of  the  safety  system.  Work  then  began  on  setting  up  the  process-control  electronics  and 
wiring. 

Combustor  Results 

Experimental  investigations  were  conducted  on  the  250  kW  natural  gas  combustor  rig,  shown 
in  Figure  2.10.  Fuel  for  the  pilot  flame  was  supplied  to  the  combustor  section  through 
the  axial  injector  on  the  bluff  body  (shown  in  Figure  2.20),  while  for  the  main  premix 
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(a)  Combustor  section  and  window  brackets. 


(b)  Guide  sleeves  and  screens. 


Figure  2.9:  Combustor  parts. 
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Figure  2.10:  Assembled  combustor  rig. 
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Figure  2.11:  Combustor  power  spectrum  under  non-reacting  156  g/s  airflow. 


combustion,  the  fuel  and  the  air  were  routed  through  a  pre-mixer  before  it  was  introduced 
in  the  combustor.  Flow  straighteners  were  provided  along  the  length  of  the  rig  to  streamline 
the  flow.  Pressure  and  temperature  measurement  ports  were  provided  along  the  length  of 
the  rig.  Preliminary  characterization  experiments  were  performed  at  atmospheric:  pressure 
and  temperature  conditions. 

For  these  experiments,  temperatures  and  dynamic  pressures  were  measured  at  50  mm 
downstream  of  the  bluff  body  face.  The  temperatTires  recorded  represent  the  qenterline  flow 
values.  Type  B  thermocouples  were  *used*for  the  flame^afea'temperatufe4  me^u?efiitiB:ts:'" 
Dynamic-pressure  measurements  were  recorded  using  differential  pressure  transducers  from 
Honeywell  Model  SX01D.  A  remote-controlled  flame-torch  was  used  to  ignite  the  pilot  flame, 
which  in  turn  ignited  the  premixed  charge.  The  preliminary  results  are  presented  here. 


Non-Reacting  Flow  The  acoustic  characteristic  of  the  combustor  in  the  presence  of 
156  g/s  (275  scfm)  airflow  is  shown  in  Figure  2.11.  As  mentioned  earlier,  the  measure¬ 
ments  were  made  50  mm  downstream  of  the  bluff  body.  The  power  spectrum  depicts  some 
pronounced  energy  at  150  Hz,  which  may  be  the  vortex  shedding  frequency. 

Figure  2.12  shows  the  acoustic-pressure  spectrum  for  the  pilot  flame  only.  A  natural 
gas  flow  of  0.15  g/s  (0.45  scfm)  was  introduced  through  the  pilot  injector  for  this  diffusion 
burning  measurement,  while  the  air-flow  was  maintained  at  156  g/s.  The  150  Hz  peak  can 
still  be  observed,  however,  a  more  significant  low-frequency  peak  appears  at  46  Hz.  The 
flame  temperature  measured  at  the  centerline  and  at  the  same  axial  location  was  1200  K 
(1700  °F).  The  corresponding  reacting  flow  image  is  given  in  Figure  2.13. 
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Figure  2.12:  Combustor  power  spectrum  for  pilot  flame  with  156  g/s  air-flow  and  0.15  g/s 
fuel-flow. 


Figure  2.13:  Reacting  flow  field  for  pilot  flame  with  156  g/s  airflow  and  0.15  g/s  fuel  flow. 


17 


2.1.  EXPERIMENTAL  ST  UDY  OF  BLUFF-BODY  COMBUSTOR ~ 


Frequency  (Hz) 


Figure  2.14:  Power  spectrum  with  combustion  at  0.4  equivalence  ratio  with  pilot  flame. 


Premixed  Charge,  <j>  =  0.4  Premix  charge  at  an  equivalence  ratio  of  <j>  =  0.4,  comprising 
148  g/s  (260  scfm)  of  air  and  3.5  g/s  (11  scfm)  of  natural  gas,  was  then  introduced  in  the  com¬ 
bustor.  The  pilot  flame  was  maintained  at  the  earlier  mentioned  values.  The  power  spectrum 
and  the  corresponding  time  signal  are  shown  in  Figure  2.14  and  Figure  2.15,  respectively. 
The  low-frequency  excitation  was  more  pronounced  with  a  slight  shift  in  frequency  to  48  Hz. 
The  flame  image  for  this  condition  is  shown  in  Figure  2.16.  The  recorded  temperature  was 
1422  K  (2100  °F). 


'  >  .  *  'J  "■ 
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Premixed  Charge,  <f>  =  0.5  The  equivalence  ratio  was  then  increased  to  <f>  —  0.5  with  air¬ 
flow  constant  and  increasing  the  fuel  flow  to  4.4  g/s  (14  scfm).  Once  again  the  pilot  flame  was 
maintained.  The  centerline  temperature  showed  an  increase  to  1478  K  (2200  °F).  The  power 
spectrum,  shown  in  Figure  2.17,  indicated  the  disappearance  of  150  Hz  oscillations,  shifting 
of  the  low-frequency  oscillation  to  38  Hz  and  the  appearance  of  sub-harmonics.  This  implied 
a  non-linear  coupling  between  the  combustor  acoustics  and  the  heat-release  oscillations.  The 
corresponding  time  signal  and  the  reacting  flow  field  are  shown  in  Figure  2.18  and  Figure  2.19. 

During  the  experiments  the  pre-mixed  flame  could  not  be  stabilized  without  the  presence 
of  the  pilot  diffusion  flame.  Efforts  are  in  progress  to  stabilize  the  flame  without  the  pilot. 
Once  that  is  achieved,  experiments  will  be  conducted  to  measure  acoustic  pressures,  tem¬ 
peratures,  heat  release  via  chemiluminescence  and  velocity  fields  via  Digital  Particle  Image 
Velocimetry  (DPIV). 
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Frequency  (Hz) 


Figure  2.17:  Power  spectrum  with  combustion  at  0.5  equivalence  ratio  with  pilot  flame. 


Figure  2.18:  Time  signal  at  0.5  equivalence  ratio  with  pilot  flame. 


Figure  2.19:  Reacting  flow  image  at  an  equivalence  ratio  of  0.5. 


2.2  Computational  Simulation  of  Combustor 


Combustor  CFD  Calculations 

The  VACCG  combustor  consists  of  a  1  in  stem  connected  to  a  bluff-body  head  (Ri  =  1  in) 
enclosed  in  a  pipe  (Ro  =  2 in).  The  head  is  configured  to  allow  for  one  axial  injector 
(R  =  0.098  in)  and  four  radial  injectors  as  shown  in  Figure  2.20.  Each  injector  is  configured 
with  an  independent  fuel-supply  and  modulation-control.  The  radial  injectors  are  spaced  at 
90°  intervals  and  located  0.5  in  upstream  of  the  bluff-body  base.  The  upstream  conditions 
are  as  follows 


Pn2  =  0.897493  kg/m3 
po2  =  0.272515  kg/m3 
PCHi  —  0.034156  kg/m3 
u  =  16  m/s 
p  =  93793.5  N/m2. 


The  above  conditions  correspond  to  an  equivalence  ratio  of  <f>  =  0.5  and  M  «  0.05. 

In  modeling  this  flow  field,  we  utilize  the  methane/air,  chemical-kinetics  model  of  Bow¬ 
man  and  Seery  [1]  Bowman  and  Seery’s  model  assumes  a  gaseous  fluid  composed  of  twelve 
species  (N2,  02,  CH 4,  C02,  CO,  H20,  H2,  H,  OH,  O,  CH3  and  HCO)  and  coupled  to¬ 
gether  by  32  elementary  reactions.  The  model  proved  to  be  more  robust  than  the  global 
mechanisms  tried  during  the  Phase  I  project.  We  continue  to  use  the  Bowman/ Seery  model 
for  the  Phase  II  calculations  and  simulate  steady  flow  through  the  combustor  both  with  and 


Figure  2.20:  Water-cooled  bluff-body  and  stem  showing  one  axial  injector  (along  the  singular 
axis)  and  two  (of  four)  radial  injectors.  Note,  the  four  larger  holes  in  the  base  region  are 
screws. 


without  axial  injection.  Chemical  reactions  are  assumed  to  advance  at  a  finite  rate.  This  is 
generally  much  less  stiff  than  assuming  equilibrium  flow,  where  reaction-rate  coefficients  are 
numerically  infinite  ( e.g .,  ~  1015). 

A  Spalart-Allmaras  [2]  one-equation  turbulence  model  is  selected  with  a  third-order, 
upwind-biased  spatial  discretization.  The  domain  is  discretized  with  two  zones  of  dimensions 


129  x  41  and  161  x  85,  respectively.  Thermodynamics  is  modeled  by  assuming  equilibrium 
translation  and  rotation  -  that  is,  the  vibrational  contribution  to  internal  energy,  which  is 
important  at  much  higher  temperatures,  is  neglected.  A  thin-layer  approximation  is  made  in 
the  viscous  formulation;  cross  derivatives  are  generally  important  for  much  lower  Reynolds- 
number  flows  and  require  a  much  higher  mesh  density  than  is  practical  for  the  duration  and 
resources  of  this  contract.  Convergence  to  a  steady  state  is  assumed  when  mass-fraction 
contours  remain  unchanged. 

Contours  of  the  steady,  pre-mixeo  simulation  without  injection  are  shown  in  Figure  2.21. 
liuiiiing  gases  are  stabiiizedrin  the  re^fCiifeMe'fr-i’Cgioii  of  the  blafT  uod^;.Tka'4einpei5ature« 


ranges  from  264.8  K  to  2623  K  and  OH  mass-fraction  is  largest  (0.99%)  near  the  base.  The 
simulation  assumes  an  axi-symmetric  centerline,  which  is  consistent  with  the  observation 
that  chemistry  reduces/eliminates  the  oscillatory  shedding  of  vortices.  Note,  however,  the 
core  flow  of  burned  gases  do  not  spread  towards  the  outer  casing  wall. 

With  pure  methane  injected  axially  at  132  m/s,  the  temperature  and  OH  mass-fraction 
contours  are  shown  in  Figure  2.22.  The  density  of  the  methane  is  pcha  =  0.648  kg/mz.  The 
peak  temperature  increases  to  2811 A  and  the  peak  OH  mass  fraction  is  0.11%.  The  axial 
injector  significantly  changes  the  OH  mass  fraction  along  the  singular  axis.  Contours  of  the 
z-component  of  velocity  and  CHz  mass  fraction  are  given  in  Figure  2.23.  Centerline  data 
and  profile  data  are  given  in  Figure  2.24. 


3-D,  Steady  Combustor  Calculation 

In  addition  to  the  two-dimensional  combustion  cases,  we  investigated  a  flow-field  simulation 
iwth  all  injectors  operating.  This  is  a  three-dimensional  flow  field  which  we  have  modeled 
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Figure  2.21:  Axi-symmetric  VACCG  combustor  without  injection.  Temperature  contours, 
streamlines  and  grid  lines  (top)  and  OH  flooded  mass-fraction  contours  (bottom). 
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Figure  2.22:  Axi-symmetric  VACCG  combustor  with  axial  injection  normal  to  bluff-body 
base.  Temperature  cons  tours,  streamlines  and  grid  lines  (top)  and  OH  flooded  mass-fraction 
contours  (bottom). 
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Figure  2.23:  Contours  of  the  x-component  of  velocity  and  CH3  mass  fraction  for  the  axial- 
injector  VACCG  combustor. 


Radial  distance,  y/R 
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(b)  Density  profiles. 


(c)  Velocity  profiles. 


(d)  Temperature  profiles. 


Figure  2.24:  Centerline  and  profile  data  for  axial-injector  combustor. 


Figure  2.25:  Temperature  field  for  the  five-injector,  three-dimensional  VACCG  combustor. 

as  symmetric  in  the  four,  azimuthal  quadrants.  Two  zones  of  sizes  129  x  41  x  33  and 
161  x  85  x  33  (593920  cells)  are  used.  Temperature  contours  and  velocity  vectors  are  shown 
in  Figure  2.25  for  the  case  of  all  five  injectors  (one  axial,  four  radial).  Methane  with  a  density 
of  0.648  kg  /m3  is  injected  through  the  five  injection  holes  at  equal  speed  of  26.46  m/s.  This 
corresponds  to  a  mass-flow  rate  'of  4121'^  10" which  is  10%  of  the  mdthane 
mass-flow  at  the  inflow.  OH  and  CHz  mass-fraction  contours  are  given  in  Figure  2.26  and 
Figure  2.27,  respectively.  Contours  of  ^-component  of  velocity  and  stream-traces  are  shown 
in  Figure  2.28.  Centerline  and  profile  data  is  shown  in  Figure  2.29.  Comparing  with  the  two- 
dimensional  simulations  (see  Figure  2.24),  combustion  occurs  much  closer  to  the  bluff-body 
base  for  the  three-dimensional  simulation.  This  is  caused  by  the  lower  injection  velocity 
(26.46  versus  132  m/s)  in  the  3-D  case.  Peak  temperatures  reach  approximately  2500  K  in 
both  cases. 

2-D,  Unsteady  Combustor  Calculation 

An  unsteady,  axi-symmetric  simulation  of  pre-mized  methane/air  flow  has  been  completed 
using  a  one-equation  Spalart-Allmaras  model.  A  steady-state  solution  has  been  used  as  an 
initial  condition  for  the  unsteady  simulation  (second-order-accurate  in  time).  The  simulation 
has  been  run  to  T  =  45ms  using  a  time  step  of  At  =  5  x  10-7  sec.  Snap-shots  of  the 
temperature  field  are  shown  in  Figure  2.30.  After  t  =  17.5ms,  the  unsteady  flow  structure 


Figure  2.26:  OH  mass-fraction  field  for  the  five-injector,  three-dimensional  VACCG  com¬ 
bustor. 
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Figure  2.27:  CH3  mass-fraction  field  for  the  five-injector,  three-dimensional  VACCG  com¬ 
bustor. 
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Figure  2.28:  Contours  of  the  x-component  of  velocity  for  the  five-injector,  three-dimensional 
VACCG  combustor. 
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fM  Density  nrofiies. 


(c)  Velocity  profiles. 


(d)  Temperature  profiles. 


Figure  2.29:  Centerline  and  profile  data  for  3-D  combustor. 
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becomes  independent  of  the  initial  condition.  The  temporal  history  beginning  at  t  =  27.5 ms 
is  shown  in  Figure  2.31.  The  temperature  distribution  at  t  =  45.0 ms  is  shown  in  Figure  2.32. 

The  velocity  and  pressure  response  is  recorded  as  a  function  of  time  at  x/R  =  1  and 
y/R  =  1  (i. e.,  one  radius  downstream  and  one  radius  above  the  centerline).  Subtracting  the 
mean  signal,  the  perturbations  are  shown  normalized  vy  the  inlet  values  in  Figure  2.33(a). 
The  mean  velocity  for  the  signal  at  is  15.28  m/s.  The  responses  of  the  mass  fractions  in 
the  chemistry  model  are  shown  in  Figure  2.33(b).  We  see  four  spikes  in  the  OH  mass- 
fraction  response  before  t  =  0.015  s  These  spikes  are  associated  with  the  transition  from  the 
initial  condition  to  the  unsteady  solution.  The  presence  of  carbon  dioxide  (CC^indicates 
combustion  and  is  necessarily  180°  out-of-phase  with  the  methane  (CH4)  mass  fraction. 

A  second  recording  is  recorded  at  x/R  =  8  to  measure  a  correlation  between  the  signals. 
The  responses  at  x/R  =  1  and  x/R  =  8  are  shown  together  in  Figure  2.34(a).  The  mean 
velocity  at  x/R  =  8  is  27.27  m/s.  A  fast  Fourier  transform  (FFT)  of  the  data  after  t  =  0.2 
is  shown  in  Figure  2.34(b).  The  FFT  reveals  the  dominant  frequency  in  the  signal  occurs  at 
158.7  Hz  with  a  sub-harmonic  at  380.9  Hz. 

Discussion  of  the  Starting  Procedure 

In  the  beginning  of  the  project,  starting  the  flame  was  computational  challenge  because  no 
hot,  ignition  source  exists  to  induce  chemical  reactions  and  subsequent  combustion.  For 
example,  in  external  aerodynamics  of  re-entry  vehicles,  chemical  reactions  of  air  species  take 
place  aft  of  a  shock  wave.  In  other  internal-flow  cases,  a  pre-heated  wall  increases  the  flow 
temperature  to  a  range  where  chemical  kinetics  become  important.  In  the  experimental 
combustor  rig,  a  spark  ignites  the  fuel/air  mixture;  however,  no  such  numerical  analog  is 
readily  apparent. 

To  circumvent  this  problem,  we  use  the  capabilities  of  QASV  to  our  advantage.  Specifi¬ 
cally,  we  can  map  a  small  sub-set  of  species  to  a  larger  model.  This  is  important  because  the 
flow  can  be  initialized  using  very  few  species  and  then  converted  to  the  full-specie  model. 
To  that  end,  we  initialize  the  flow  with  a  cold-flow  air/methane  mixture  on  a  coarse  mesh. 
This  is  a  very  rapid  process  which  establishes  a  good  initial  guess  to  the  cold-flow  simulation 
-  /or  a  small  numb  or-  of  -points-  oz-i  species.  Needless.to-.^jayy^ie:G@m|5.'ata;tiorial.;> 

time  required  for  a  full-specie  model  on  a  fine  mesh  would  be  substantially  larger. 

With  the  density,  velocity  and  pressure  fields  well-defined  for  the  correct  equivalence 
ratio  and  geometry,  we  perform  two  conversions.  First,  we  convert  the  air  density  to  79% 
nitrogen  and  21%  oxygen  (by  mole),  while  copying  the  methane  density  unchanged.  There 
are  12  species  in  the  Bowmann/Seery  model,  and  the  remaining  specie  densities  in  the  model 
are  set  to  zero.  Secondly,  the  temperature  in  the  combustion  zone  is  increased  artificially  by 
lowering  the  specie  densities  while  holding  the  pressure  constant.  Because  T  =  p/ (^  PiRi), 
a  decrease  in  all  the  densities,  pi,  results  in  an  increase  in  the  temperature,  T.  We  decrease 
the  densities  by  a  factor  of  ten  which  gives  a  combustion-zone  temperature  in  the  range  of 
the  adiabatic  flame  temperature.  This  solutin  is  then  used  with  chemical  kinetics  active  until 
the  flow  reaches  a  physically  realistic  solution. 

The  experience  with  numerically  “igniting”  internal  combustor  flows  has  demonstrated 
the  need  to  place  restrictions  on  the  magnitude  of  the  chemical  source  terms.  Many  methods 
for  starting  the  chemistry  were  tried,  but  only  the  above  procedure  proved  to  be  well  behaved. 
Often  times,  the  simulation  quickly  diverged  as  the  mixture  progressed  to  unrealistic  states 
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(a)  t=27.5ras 


(b)  t =30.  Oms 


(c)  t=32.5ms 


(d)  t=35.0ms 


(e)  t=37.5ms 


(f)  t=40.0ms 


(g)  t=42.5ms  (b)  t=45.0ms 

Figure  2.31:  Temperature  contours  for  axi-symmetric,  two-dimensional  blunt-body  combus¬ 
tor  without  methane  injection.  Sequence  shows  temporal  independence  from  initial  condi¬ 
tion. 
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Normalized  x  coordinate,  X/R 


Figure  2.32:  Temperature  profile  at  y/R  =  1. 


( e.g infinite  temperatures,  negative  densities  or  pressures).  One  procedure  for  igniting 
the  flow  that  was  very  stiff  involved  increasing  the  pressure  for  a  small  number  of  cells 
in  the  recirculation  region.  This  resulted  in  an  increase  in  the  temperature,  but  caused  a 
tremendous  disturbance  in  the  velocity  and  pressure  fields.  As  a  result,  the  maximum,  stable 
time  step  was  a  very  small  value  on  the  order  of  A t  =  10~9  sec.  This  necessity  brought  about 
the  invention  of  a  procedure  that  we  call  “Damkohler  limiting” .  As  a  result  of  this  work,  an 
AIAA  paper  was  produced  which  will  be  presented  at  the  43rd  Aerospace  Science  Meeting 
on  January,  2005.  We  include  below  an  excerpt  from  that  paper  in  this  report. 


Nfcss  Fraction  at  x/R=l,  y/R=l 
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Time  (sec) 

(b)  Mass  fraction  response. 


Figure  2.33:  Time  history  of  ^-component  of  velocity  and  mass-fractions  at  x/R  =  1  and 
y/R  =  1. 
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Frequency  (Hz) 

(b)  FFT  of  signal. 

Figure  2.34:  Time  history  of  ^-component  of  velocity  at  x/R  =  1  and  x/R  =  8  and  y/R  =  1. 
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(Excerpt  from  AIAA  Paper  2005-1400) 


The  Damkohler  Number 


The  Damkohler  number  is  generically  defined  as  the  ratio  of  the  fluid-dynamic  time  scale  to 
a  second  physically  relevant  time  scale  such  as  chemical  kinetics,  binary  diffusion  or  heat  and 
momentum  transfer.  For  example,  we  can  consider  the  inverse  of  the  cell  Reynolds  number 
as  a  Damkohler  number  for  viscous  diffusion  as  follows 


Da^  —  1  j  Re  Ax  — 


a  (MSI 

pUAx  (pAx2/p.) 


(2.1) 


For  a  viscous  flow,  diffusion  of  momentum  becomes  important  when  the  viscous  Damkohler 
number  ( i.e .,  the  inverse  of  the  Reynolds  number)  increases.  More  conventionally,  diffusion 
plays  a  more  important  role  as  the  Reynolds  number  decreases.  A  large  Damkohler  number 
means  that  the  time  scale  of  interest  is  much  smaller  than  the  fluid-dynamic  time  scale.  In 
other  words,  a  small  time  scale  means  that  the  physical  process  associated  with  the  scale  is 
happening  quickly  compared  to  the  convection  of  fluid.  Likewise,  a  large  time  scale  means 
the  process  is  progressing  slowly.  This  has  particular  relevance  for  chemically  reacting  flows 
where  the  Damkohler  number  is  defined  as 

Daj,  =  Hi.  (2.2) 

T  i 


A  large  Damkohler  number  (Da  »  1)  means  that  the  chemical  time  scale  is  much  larger 
than  the  fluid-dynamic  scale  and  the  flow  is  tending  toward  equilibrium.  Likewise,  a  small 
Damkohler  number  means  the  flow  is  nearly  frozen  since  the  chemical  changes  would  be 
lagging  far  behind  the  fluid-dynamic  changes.  The  local  fluid-dynamic  time  scale,  7/^  is 
defined  here  as  the  time  required  for  a  fluid  particle  to  advect  across  one  cell. 


Numerical  Stiffness 

When  a  physical  system  has  drastically  different  characteristic  times,  numerical  stiffness 
results  which  limits  the  progression  to  a  steady  state  by  the  smallest  time,  ij.s-.seen  in  the 
s^ar,  adveetioii  equation,  the  affect  of The  source  -term7 constrain  the  time  step  to 
min(aAx,  2 r).  In  finite-rate  chemical  kinetics,  rch  — ►  0  as  a  mixture  approaches  equilibrium 
for  a  fixed  flow  speed.  Therefore,  rapid  chemical  kinetics  causes  numerical  stiffness  because 
the  numerical  time  step  must  approach  zero  in  order  to  maintain  stability.  This  is  undesirable 
for  steady-state  problems  and  Godfrey  [3]  shows  that  implicit  algorithms  bypass  this  transient 
stiffness  by  re-scaling  the  chemical  characteristic  time. 

Eberhardt  and  Imlay  [4]  show  thaf^the  fastest  reaction  slows  all  reactions,  and  thus  slows 
down  steady-state  convergence.  Any  time  scale  that  is  shorter  than  the  fluid-dynamic  scale 
is  considered  a  “sub-grid”  time  scale.  The  fluid-dynamic  scale  is  the  time  required  for  a 
fluid  particle  to  advect  across  a  cell.  If  the  chemical  composition  is  changing  so  fast  that  the 
mixture  reaches  equilibrium  before  crossing  a  cell,  then  the  reaction  rates  can  be  reduced  to 
delay  the  rate  at  which  equilibrium  is  reached.  In  other  words,  there  is  no  need  to  resolve 
chemical  time  scales  which  are  smaller  than  the  grid  can  resolve.  Eberhardt  and  Imlay  use 
a  maximum  Damkohler  number  of  ten  (Damax  =  10)  and  remark  that  the  effect  is  to  smear 
shocks  slightly.  We  will  see  in  our  results  that  the  maximum  error  introduced  by  the  limiting 
occurs  at  a  shock  discontinuity. 
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Chemical  Time  Scales 

For  unsteady  flow  computations,  the  Damkohler  number  for  chemistry  is  defined  as  the  ratio 
of  fluid-dynamic  to  chemical  time  scales.  For  the  remainder  of  this  paper,  the  chemical 
Damkohler  number  will  be  denoted  as 


Da  _  _min{r/d,  At} 


Tch 


(2.3) 


Generally,  the  chemical  time  scale  is  inversely  proportional  to  the  fastest  reaction  rate  in 
the  model.  Damkohler  limiting  modifies  the  fastest  reaction  rate  while  maintaining  the 
equilibrium  constant  (i.e.,  the  backward  rate  is  changed  to  compensate  for  the  change  in  the 
forward  rate  so  that  their  ratio  remains  the  same). 

We  can  represent  a  general  chemical-kinetic  system  composed  of  N  species  and  J  reac¬ 
tions  as  follows 

v'ljXi  +  v'i'jXi  - V  v’njXn  #  v"l,jX\  +  v" 2,3X2  H - b  1 '"njXn 

j  =  1,2,...,  J,  (2.4) 

where  1/  ■  and  v"  -  are  the  stoichiometric  coefficients  of  species  X{  in  the  jth  reaction.  The 

ijJ  li3 

rate  of  production  of  species  s  is  given  as 


dp 

dt 


N 


N 


1  -  Mi  SKi  -  Kj)  [kf<j  kb'j  ft (Xft  ^‘,3 

j—\  1= 1  1=1 


i  =  1,2,  ...,N ,  (2.5) 


where  kfj  and  kbj  are  the  forward  and  backward  reaction  rates.  The  forward  rates  are 
determined  from  the  Arrhenius  equation 


kf  =  CTV  e-£/fcr, 


(2.6) 


and  the  backward  rates  are  determined  from  the  ratio  of  the  forward  rate  to  the  equilibrium 
constant 


-  -  kb'  ■ 


kf 

Kp 


Jo  7 


7) 


The  equilibrium  constants  are  determined  from  either  the  Arrhenius  equation,  equilibrium 
curve  fits,  or  from  the  McBride  Curve  fits  and  the  minimization  of  Gibbs  free  energy. 

Bimolecular  Reaction  We  are  interested  in  the  chemical  time  scale  for  the  elementary 
reactions  that  compose  a  chemical  kinetics  model.  Following  Turns  [5],  we  consider  the 
following  bimolecular  reaction 

A  +  B  —>  products  (2.8) 

and  its  rate  equation 

M.  =  -k[A][B].  (2.9) 

at 

We  substitute  the  defect  variable,  x  =  [A]  —  [A]o  =  [B ]  —  [B]o,  into  the  rate  equation  and 
integrate  to  yield 

dx  =  -kt  +  C. 


/ 


(x+  [A]o)(x+  [B]o) 


(2.10) 
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Reaction 

Example 

Tch 

2  A  — ►  products 

20  — >  O2 

e  —  1 

2[A]0  k 

A  +  M  — >  products  +  M 

H20  +  M-^H  +  OH  +  M 

1 

[M]o  k 

2  A  — ►  products  +  A 

2  O2  — ►  2  0  +  O2 

e  —  1 
[A]0k  • 

A  +  B  —*  products 

H2  +  02->2H0 

C(x) 

(Wo  -  [Bio)  k 

A  +  B  +  M  — *  products  +  M 

N  +  O  +  M^NO  +  M 

C(x) 

(Wo  -  |B]o)  k\M\ 

Table  2.2:  Chemical  time  constants  where  x  =  [A]o/[jB]o  and  [A]  is  assigned  to  the  smaller 
concentration. 


The  left-hand  integral  can  be  determined  to  be 

_ _ = _ I _ in  (W 

(®  +  [i4]0)(*  +  [B]o)  [A]o-[B}0  V[A] 

and  setting  this  result  equal  to  the  right-hand  side  of  Eqn.  (2.10),  we  can  solve  for  the 
temporal  decay  of  [A]/[B]  as  follows 

jjj  =  jgjj;  <*p{(Wo  -  [B]„) «}.  (2.12) 

This  expression  is  used  to  determine  tficrfime  of  the  ’spdcie  coiicetf- ' 

trations. 


(2.11) 


Definition:  Time  Constant  The  time  constant  is  defined  as  the  time  elapsed  for  an 
initial  concentration  to  decay  by  1/e.  Assigning  [A]o  to  be  the  smaller  of  the  two  concen¬ 
trations,  we  solve  Eqn.  (2.12)  for  tch  by  solving  for  the  time  when  [A]/[A]o  =  1/e  where  we 
use  the  fact  that  [B\/[B\o  =  ([A]0/[B]0)([A]/[A]o  -  1)  +  1.  We  obtain  the  time  constant  for 


a  bimolecular  reaction  as  . 

£{[A]0/[f?]o} 

“*  ([A]o  -  |B]o)  * 


(2.13) 


where 


C{[A}o/[B}o}  =  In 


[A]/[A]  o 


[B]/[B]  o 


=  -  In 


[A]/[A]o=l/e 


(2.14) 


Table  2.2  lists  the  chemical  time  scales  for  several  common  elementary  reactions. 
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A  General  Method  for  Computing  the  Chemical  Time  Constant  As  demonstrated 
above,  the  analytical  expression  for  the  chemical  time  constant  depends  on  the  form  of  the 
elementary  reaction.  Rather  than  devise  expressions  for  every  possible  reaction,  we  seek  a 
general  method  which  applies  to  every  case.  In  the  context  of  Eqn.  (2.4),  consider  the  jth 
chemical  reaction  in  a  general  chemical-kinetics  model 

v'i,jX i  +  v' 2,3X2  H - 1-  v'njX n  #  v"ijX\  +  1 '"2,3X2  -\ - h  v"njXn-  (2.15) 


The  defect  law  for  each  specie  can  be  written  as  follows 


[Xi]  -  [Xilo  [X2]  -  [Agjo 


vh  -  vi,i 


[XN]  ~  [Xn]o 
vn,o-u'n,o 


(2.16) 


We  can  compute  the  chemical  time  scale  numerically  by  solving  an  ordinary  differential 
equation  for  the  specie  with  the  smallest  negative  defect.  That  is,  we  assume  that  v\yj  = 
v,{j  —  v[j  <  0  and  that  the  species  are  listed  in  order  of  increasing  defect 

pik<Ho<...<Eid£.  (2,i7) 

Kjl  K  j  |  Wnj\ 

Under  these  conditions,  we  can  solve  the  following  ODE  for  the  decay  of  [Xi] 


d\X 1] 

dt 


=  vld  kfjlXifa [Xjfis  ■  •  •  \Xitf"'*- 


(2.18) 


The  defect  law  (Eqn.  (2.16))  can  be  used  to  write  all  the  specie  concentrations  if  Eqn.  (2.18) 
in  terms  of  [Xi].  A  similar  ODE  can  be  written  for  the  smallest  decreasing  specie  in  the 
backward  reaction  with  rate  coefficient  . 

To  solve  Eqn.  (2.18)  for  the  chemical  time  scale  we  use  a  second  order  mid-point  method 
to  determine  the  time  when  the  specie  decays  by  a  ratio  of  1/e.  We  first  make  a  first-order 
approximation  using  the  Euler-explicit  method  for  an  initial  guess 


[*!]*  =  [Xl]0  + 


d[Xx) 


dt 


*-* 

^ ch" 


(2.19) 


Solving  for  r*h  and  noting  that  [Xi]*/[Xi]o  =  1/e,  we  have 


T*h  = 


(1/e  -  1). 


(2.20) 


The  second-order  approximation  can  be  written  by  re-evaluating  the  right  hand  side  of 
Eqn.  (2.18)  using  the  mid-point  concentration  (i.e.,  [Xi]o  +  [Xi]q  (1/2  r*h)).  Our  chemical 
time  scale  for  any  elementary  reaction  can  then  be  written  as 
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Fluid  Dynamic  Time  Scales 

The  fluid-dynamic  time  scale  is  determined  locally  for  each  cell  in  the  flow  field.  For  a 
structured  mesh  in  a  finite-volume  setting,  we  compute  a  time  scale  for  each  coordinate 
direction  and  choose  the  smallest  one  in  magnitude.  For  our  calculations,  we  compute  the 
following  three  time  scales 

_ 2  A  Vijk _ 

l^i— 1/2 1  Aj4,— 1/2  +  l^i+l/2|AAi+1/2 

2  AVjjfe _ 

l“j— l/2|AAj_l/2  +  I^j+1/2|Aj4j+1/2 

_ 2  A  Vijk  _ 

l^fc— I/2I  AAfc_!/2  +  |%fl/2|A^A:+l/2 

where  u  =  V  •  n  is  the  velocity  component  normal  to  each  cell  face.  The  fluid-dynamic  time 
scale  is  the  minimum  of  the  three  possibilities. 


Tfdi  = 
Tfdj  ~ 

Tfdk  = 


Application 

For  each  reaction  (j)  in  a  model,  we  determine  Damkohler  numbers  based  on  the  forward  and 
backward  time  constants.  If  the  larger  Damkohler  number  exceeds  a  threshold  (i.e.,  Damax), 
then  the  corresponding  reaction  rate  is  modified  in  such  a  way  as  to  maintain  the  equilibrium 
constant.  That  is, 

•  If  Da jj  >  Dabj  and  Da fj  >  Damax,  then  the  modified  rate  coefficient  becomes 


minlr/d,  At} 


(2.22) 


•  If  Da^j  >  Dafj  and  Da^j  >  Damax ,  then 


u 


min  {Tfd,  At}) 


If  limiting  is  applied  to  either  rate  coefficient,  the  opposite  rate  coefficient  is  adjusted 
maintain  the  following  equality 


Because  the  method  maintains  the  original  law  of  mass  action,  the  equilibrium  distri¬ 
bution  remains  unchanged  and  only  the  rate  to  reach  equilibrium  altered. 


In  the  above,  u  =  rchk ,  which  makes  u  dependent  on  the  type  of  reaction.  However,  we  can 
assume  that  bimolecular  reactions  will  occur  faster  than  termolecular  ones.  If  we  assume 
that  one  concentration  in  the  bimolecular  reaction  is  much  larger  than  the  other,  then  u> 
becomes  in  the  limit 

^  max{[A]0,[B]o}’ 


(2.25) 
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Supersonic  Air  Nozzle 

We  consider  steady,  supersonic  flow  of  a  five-species,  expanding  air  mixture  through  a  nozzle 
with  the  following  inlet  conditions: 

pNi  =  0.00048557  kg/m3 
po2  =  1.55682  x  10'7  kg/m3 
pNO  =  9.9556  x  10~6  kg/m 3 
pu  =  0.0146888  kg/rn3 
po  =  0.00460352  kg/m3 
U  =  4000  m/s 
T  =  9000  K. 

At  this  high  temperature,  the  inlet  gas  is  composed  of  primarily  dissociated  nitrogen  ( pn / P  = 
74.2%)  and  oxygen  (p0/p  =  23.2%)  with  a  small  amount  of  diatomic  nitrogen  (2.4%)  and 
trace  amounts  of  O2  and  NO. 

We  simulate  the  nozzle  flow  using  chemical-kinetics  model  of  Kang  et  al.  [6]  on  a  11  x  41 
mesh  with  uniform  spacing.  Time  integration  consists  of  a  first-order  implicit  marching 
scheme  until  the  relative  residual  tolerance  declines  five  orders  of  magnitude.  The  distribu¬ 
tion  of  Damkohler  limiting  is  shown  in  Figure  2.35(a).  The  contours  show  that  the  Damkohler 
number  is  largest  at  the  in-flow  plane  (Da  =  761.5)  where  the  temperature  (and  reaction 
rates)  are  largest  and  the  flow  speed  is  the  smallest.  As  the  flow  expands,  the  tempera¬ 
ture  decreases  (see  Figure  2.35(b))  and  the  flow  speed  increases.  Therefore,  the  Damkohler 
number  decreases  downstream  until  Da  =  42.7  at  the  exit. 

In  order  to  examine  the  effect  of  Damkohler  limiting  on  solution  quality,  a  number  of 
Damkohler-limited  cases  are  run  with  Damax  ranging  from  Damax  :  0.1  —*  10.  In  Figure  2.36, 
we  present  the  mass-fraction  distributions  for  the  various  Damkohler  numbers  compared  to 
the  unlimited  case.  For  all  Damkohler  numbers  except  Damax  =  0.1,  the  results  are  virtually 
identical  between  the  unlimited  and  limited  cases.  In  the  most  severely  limited  case,  the 
primary  source  of  the  error  occurs  at  the  inlet.  The  trend  in  the  mass  fraction  downstream 
of  the  inlet  closely' matches'  the' unliiriited'?caser  We" Emphasize  that'  'the "fh'ost  severemase’  - 
limits  the  reaction  rates  by  a  factor  of  7615  smaller  than  the  original  model.  However,  the 
solution  remains  virtually  unchanged.  While  the  errors  in  Figure  2.36(b)  appear  to  be  large, 
the  scale  of  the  axis  is  very  small  and  the  peak  error  is  less  that  0.05%  at  Damax  =  1.0. 

(End  of  Excerpt  from  AIAA  Paper  2005-1400) 

2.3  Implementation  of  Sensitivity-Equation  Method  in  GASP 

Inviscid  Sensitivity  Modifications 

One  objective  of  the  Phase  II  portion  of  this  project  is  to  incorporate  the  sensitivity  software 
that  has  been  developed  in  the  stand-alone  sensitivity  solver,  SSAfSS,  into  AeroSoft’s  com¬ 
mercial  flow  solver,  QASV-  The  reason  for  this  transition  is  that  time-accurate  calculations 
require  a  tight  coupling  between  the  flow  and  sensitivity  solvers.  The  stand-alone  capability 
(which  is  so  beneficial  for  steady-state  calculations)  is  still  available  through  QASV,  but 
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Figure  2.37:  GASP  graphical  user  interface  displaying  the  diamond-airfoil  sensitivity  solu¬ 
tion. 

many  features  that  have  been  developed  for  the  flow  solver  can  also  be  applied  to  sensitivity 
problems.  These  include  the  user  interface,  distributed  parallel,  overset  methods,  implicit 
and  explicit  temporally  accurate  methods,  and  many  thermodynamic  and  chemical  reaction 
models. 

As  a  way  of  beginning  the  testing  under  the  QASV  framework,  an  inviscid  flow  and 
sensitivity  calculation  were  conducted  for  supersonic  flow  over  f  diamond  airfoil.  The  user 
interface  for  this  problem  is  shown  in  Figure  2 ,37.  where,  the  design  variable  is  the  fre^stream 
velocity.  From  the  interface,  the  user  can  output  data  of  four  different  types:  flow  data, 
sensitivity  data,  near-by  data,  or  uncertainty  data.  Very  few  modifications  we  required  by 
the  user  to  transform  the  GUI  from  running  a  flow  problem  to  running  a  sensitivity  problem. 

The  pressure  distribution  along  the  airfoil  for  a  5%  perturbation  from  Moo  =  2.0  to 
Mm  =  2.1  is  shown  in  Figure  2.38.  The  near-by  approximation  compares  well  with  a  new 
calculation  at  Mm  =  2.1. 

Viscous  and  Turbulent  Sensitivity  Modifications 

A  subsonic,  flat-plate  boundary  layer  flow  is  used  to  verify  that  the  sensitivity  equation 
method  is  being  utilized  correctly  for  the  eddy-viscosity  models  of  QASV.  Each  of  these  mod¬ 
els  has  been  tested  and  incorporated  into  AeroSoft’s  stand-alone  sensitivity  solver,  SENSE. 
However,  because  the  discretization  methods  in  QASV  are  slightly  different,  the  coding  task 
was  not  as  straight-forward  as  importing  the  pre-existing  code.  After  completing  coding 
for  viscous  fluxes  and  turbulence  source  terms,  the  sensitivity  results  for  three  turbulence 
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Diamond  Airfoil  Sensitivity  to  Free-stream  Velocity 

Near-by  Flow  Approximation  at  M=2.1 
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Figure  2.38:  Near-by  pressure  approximation  compared  to  a  new  calculation  to  for  the 
diamond  airfoil  sensitivity  to  free-stream  velocity.. 


models  are  shown  in  Figure  2.39.  These  results  are  consistent  with  what  has  been  obtained 
using  SENSE  in  Godfrey  and  Cliff  [7]. 


Time- Accurate  Numerics 

The  final  modification  needed  to  implement  the  sensitivity  software  in  QASV  involved  the 
time-marching  schemes  Rensitiyity.madificatione 

use  the  explicit  and  implicit  time-accurate  algorithms.  This  process  required  changing  data 
structures  from  manipulating  only  flow-related  data  to  manipulating  general  sets  of  data 
(flow  and  sensitivity).  As  such,  a  level  of  abstraction  had  to  be  added  to  the  solution-related 
data  structures.  Being  written  in  C++,  QASV  is  easily  adapted  to  such  an  abstraction.  With 
the  time-accurate  capability  complete,  a  user  can  run  time-accurate  sensitivity  calculations 
in  conjunction  with  either  fixed  or  time-accurate  flow  solutions. 

As  a  verification  case,  the  variable  accuracy  Runge-Kutta  and  quasi-Newton  algorithms 
have  been  applied  to  a  typical  shock-tube  case  where  the  design  variable  is  the  high-pressure- 
region’s  pressure  (see  also  Appel  [8]).  The  sensitivity  results  for  various  time-integration 
schemes  are  shown  in  Figure  2.40  with  comparison  to  shock-tube  theory.  This  problem  was 
run  by  ICAM  during  the  Phase  I  project  and  is  now  being  reproduced  for  the  first  time  with 
the  newly  developed  QASV  capability.  Since  the  results  are  in  very  good  agreement  with 
the  theoretical  sensitivities  (except  at  the  discontinuities),  we  are  confident  that  the  schemes 
have  been  implemented  correctly  for  time-accurate  sensitivity  calculations.  Thus,  we  could 
move  forward  toward  investigating  the  combustion-related  sensitivity  problems. 
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Turbulent  Flat  Plate  Sensitivity  to  Edge  Velocity 

ReL=2xl06,  M_=0.3 


Figure  2.39:  Sensitivity  results  for  a  turbulent  flat  plate. 

2.4  Feedback  Control  of  1-D  Model  Problem 

This  section  of  the  report  is  similar  in  style  to  a  technical  paper  rather  than  a  step-by-step 
summary  of  the  activities  carried  out.  We  describe  here  the  work  completed  in  controlling  a 
model,  acoustic-instability  problem. 

One-dimensional  Non-linear  Flow  Equations 

In  Venkateswaran  et  al.  [9]  and  Wang  et  al.  [10],  the  authors  use  a  boundary  condition  which 
couples  the  pressure  and  velocity  perturbations  to  cause  an  initial  acoustic  wave  to  become 
unstable.  A  one-dimensional,  pressure  perturbation  grows,  at  first  linearly,  until  reaching  a 
non-linear  limit  cycle.  We  wish  fcx  reproduce, th^e  results  and  then  consider  some  control 
methods  to  suppress  the  instability- 

For  generality,  consider  the  one-dimensional,  unsteady  flow  of  a  compressible,  perfect 
gas  in  a  tube  of  constant  area.  Conservation  of  mass,  momentum  and  energy  lead  to  the 
following  partial  differential  equation  (PDE) 


+  dF(Q(t>_g))  =  s(t)X)|  o  <  a:  <  L, 
at  ax 


t>  0, 


(2.26) 


where 


Q  = 

In  these  equations 

•  p(t,  x)  is  the  density, 

•  u(t,  x)  is  the  velocity, 


p 

pu 

0 

pu 

,  F  = 

pu2  +  p 

,  and  S  = 

0 

.  m . 

puho 

_  Se(t,x)  _ 
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Sod’s  Problem 

10:1  Pressure,  8:1  Density,  T=10’3,  At=10^,  ENO  limiting 


-1  -0,8  -0.6  -0.4  -0.2  0  0.2  0.4  0.6  0.8  1 


X 


(a)  Density  sensitivity. 


Sod’s  Problem 

10:1  Pressure,  8:1  Density,  T=10 3,  At=10^,  ENO  Limiting 


x 


(b)  Pressure  sensitivity. 


Figure  2.40:  Sensitivity  solution  to  the  10:1  shock  tube. 
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Figure  2.41:  Initial  triangular  pressure  pulse  for  one-dimensional  control. 


•  e(t,  x )  is  the  internal  energy  per  unit  mass, 

2 

•  eo(t,  x)  is  the  total  energy  per  unit  mass  (=e  +  \), 

•  p(t,x)  is  the  pressure  (p  =  (7  —  l)pe), 

•  ho(t,x)  is  the  stagnation  enthalpy  per  unit  mass  (=  eo  +  *■),  and 

•  Se(t,x )  represents  input  of  thermal  energy  (e.g.  watts/m3). 


Here  7  is  the  usual  ratio  of  specific  heats. 

We  are  interested  in  the  growth/control  of  a  triangular,  pressure  pulse  shown  in  Fig¬ 
ure  2.41.  The  initial  flow  state  consists  of  a  uniform,  subsonic  flow  with  a  ±1%  perturbation 
in  the  pressure  at  the  mid-section  of  the  pipe.  Tie  length  and  conditions  of  the  pipe  are 

consistent  with  those  e^merieneed  int>li6iiACCG  comb  v,ryor....;c.. w  ^ 


L  =  1.4732  m 
po  —  1.2042  kg/m3 
uo  =  15  m/s 
po  =  93760  N/m2 


70  =  1.396 


The  inflow  boundary  condition  described  earlier  couples  pressure  and  velocity  perturba¬ 
tions  to  sustain  an  aero-acoustic  instability.  The  boundary  condition  can  be  conceptualized 
as  a  collapsed-combustion  zone  where  pressure  perturbations  lead  to  in-phase  velocity  per¬ 
turbations.  This  boundary  condition  is  a  simplification  of  the  role  played  by  combustion 
in  the  instability  mechanism.  The  coupled  boundary  condition  [10]  satisfies  the  following 
relationship: 


1  + 


•e 


PQUQ 

Pi 


-  Uq 


(2.27) 
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where  we  have  selected  a  =  1  for  all  the  calculations.  'The  state  at  the  boundary  face  is 
given  the  subscript  “1”.  The  “0”  subscript  denotes  a  reference  state  which  is  the  constant 
flow  about  which  we  make  the  linearization.  Perturbations  are  denoted  with  a  single  prime, 
( e.g v!  =  u  -  uQ).  To  complete  the  boundary  condition,  pressure  disturbances  are  extrapo¬ 
lated  from  the  interior  and  the  entropy  (p/p7)  is  held  constant.  The  outflow  boundary  serves 
as  a  node  for  the  pressure  wave.  All  variables  are  extrapolated  except  the  pressure,  which  is 
held  constant. 


Linear  System 

Assuming  smooth  solutions,  the  quasi-linear  form  of  Eqn.  (2.26)  is  written 

9Q(t,x)  argq  = 

at  dQdx  K  ’ 


where 


0  1  0 
-(3~7)t  (3~7)_u  7-1 

[  (7  —  l)u3  —  que 0  qeo  —  33ylti2  7 u  j 

We  assume  an  underlying  uniform,  steady  flow  and  write 


dF 

dQ 


(2.28) 


Q(t,x)  —  Qo  ■+■  Q  (t,^), 

where  Qo  is  the  constant  vector  value  of  the  conserved  quantities  in  the  unperturbed  flow, 
and  Q(*,:r)  is  a  perturbation.  The  temporal  history  of  Q (t,x)  is  described  by  the  linear 
PDE 

(2.29) 


dQ(t,  x)  dF 

at  +  dQ 


§  =  S(t,x). 


Primitive  Variable  Formulation 


For  several  purposes,  such  as  the  discussion  of  boundary  conditions,  it  is  of  interest  to 
consider  s.  change  .of  dependent  variable&from  the  consor-vative  set  (Q)  to  so-called  primitive 
variables.  Here  we  choose 

r  p  1 


where,  as  noted  above,  p  is  the  density,  u  is  the  velocity,  and  p  is  the  pressure.  It  can  be 
readily  verified  that  the  conservative  and  primitive  variables  are  related  according  to: 


Qi 

qi 

q(Q)  = 

Q2/Q1 

and  Q(q)  = 

qiq2 

.(7-1)  [Q3  -  Qi/(2Qi)j  _ 

.^17  +  jqiql. 

The  Jacobians  of  these  transformations  are  given  by 


dQ 

dq 


1  0  0 

up  0 

u2 /2  pu  I/7  —  1 


t 


(2.30) 
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and  by 

1  0  0 
-u/p  1/p  0 

(7  -  l)tt2/2  -(7  -  l)tx  7  -  1 

The  linearized  PDE  for  the  primitive  variables  is 


dq 

dQ 


(2.31) 


dqjt,  x )  dq 

at  dx 


s  (t,x), 


(2.32) 


where 


u  p  0 

0 

A  = 

0  u  1/p 

,  s  = 

0 

0  7 p  u 

.  (7-  l)Se{t,x)  . 

(2.33) 


Notice  that  the  linearized  system  in  primitive  variables  take  the  form  of  a  state-space 
system 

w  =  [A]w  +  [B]u  (2.34) 


where  w  represents  perturbations  in  the  primitive  flow  variables  and  u  is  a  vector  containing 
zeros  and  the  control  variable,  which  we  will  discuss  later. 


Uncontrolled  Response 

In  the  uncontrolled  case  ( i.e Se  =  0),  the  temporal  response  of  the  pressure  at  the  inflow 
boundary  is  shown  in  Figure  2.42(a)  for  the  both  the  non-linear  (Eqn.  (2.26))  and  linear 
(Eqn.  (2.32))  systems.  The  envelope  of  the  response  is  given  in  Figure  2.42(b)  which  shows 
the  peak  values  of  each  oscillation.  For  the  remainder  of  the  report,  we  use  the  envelopes  of 
the  response  for  clarity  and  to  reduce  the  amount  of  ink  in  the  figures. 

The  linear  and  non-linear  formulations  agree  well  with  one  another  until  the  limit  cycle 
begins.  The  linear  system  predicts  continued  growth  whereas  the  non-linear  terms  cause  a 
limit  cycle  at  t  «  1.6  sec.  The  magnitude  of  the  limit  cycle  is  approximately  p/po  =  6.3% 
with  a  frequency  of  54.93  Hz  (approximately  4L/ao). 

Control  Method 

We  postulate  a  control  in  our  problem  which  represents  the  input  of  thermal  energy  on  the 
interior  of  our  spatial  interval  [0,  L] .  The  mathematical  model  is  represented  as 

Se(t,  x)  =  vh(t)bh(x),  (2.35) 

where  b h  is  a  fixed  spatial  function,  while  the  control  amplitude  is  modulated  by  the  time 
function  Vh .  Our  main  purpose  is  to  study  the  possibility  of  stabilizing  the  system  through 
appropriate  feedback  control. 

The  bh  term  represents  heat  addition  and  is  loosely  associated  with  injected  methane  fuel 
in  the  combustor  problem.  The  heat  addition  is  distributed  in  space  according  to  a  standard 
normal  distribution  the  Gaussian  distribution)  which  can  be  written  as 

1  £2 

= -^=  exp{— ^ },  where  £  = 


a 


(2.36) 


Pressure  fluctuation,  pvp _  Pressure  fluctuation,  p Vp^ 
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Pressure  Response  to  a  1-D  Bombing  Pulse,  p’=0 

&p/P'=[%,y=1.135I,  M=0.3, 2nd-Order  RK, 


Time  (sec) 

(a)  Linear  and  non-linear  response. 


VACCG  Pressure  Response  Envelope  to  a  1-D  Bombing  Pulse 

Ap/Pt=\%,y=L396,  A/=0.05, 2nd-Order  RK,  A  t=10‘5 


(b)  Envelopes  of  response. 


Figure  2.42:  Pressure  response  and  envelope  for  both  the  linear  and  non-linear  systems. 
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Figure  2.43:  Gaussian  distribution  used  in  model  for  combustion  control. 

This  distribution  is  shown  in  normalized  form  in  Figure  2.43.  The  mean  of  the  distribution 
is  represented  by  /i  and  the  variance  by  a;  the  integral  under  the  Gaussian  distribution  is 
unity.  Our  control  is  chosen  such  that  the  total  amount  of  heat  added  to  the  system  is  a 
fraction  of  the  amount  necessary  to  produce  thermal  choking.  That  is, 

bh(O  =  eAh*0G&H,a).  (2.37) 

In  Eqn.  (2.37),  A hi  =  h^-ho  is  the  amount  of  heat-addition  required  to  produce  thermally 
choked  flow.  From  Anderson  [11],  the  amount  of  heat  addition  needed  for  choking  can  be 
determined  from 

h0  2  (7  +  1  )M2 
hi  ~  (1  • 7M2)2 

We  choose  e  =  10-3  (i.e.,  0.1%  of  choked 'heaTad'ditlbh)  and  the  control  is  applied  at  the 
front  of  the  domain  by  selecting  /. 1  =  L/ 10  and  cr  =  L/ 60. 

Abstract  LQR  control  problem 

Our  objective  is  to  stabilize  the  non-linear  system  described  by  the  Euler  equations  subject  to 
the  coupled  boundary  conditions.  One  way  to  construct  a  stabilizing  controller  is  to  formulate 
an  optimal-control  problem  with  linearized  dynamics  and  a  quadratic  cost  functional.  The 
dynamics  are  linearized  about  a  uniform  steady  flow  and  our  cost  functional  is  defined  using 
inner  products  as 

J(u)  =  J  |{(/  wtQw dx^j  +  utRu|  dt,  (2.39) 

where  we  select  Q  =  diag {Qp,  Qu-,  Qp}  to  be  a  diagonal  matrix  of  prescribed  non-negative 
functions  that  penalize  perturbations  in  density,  velocity,  and  pressure,  respectively.  R  is 


1  +  1—±M2 


(2.38) 
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a  matrix  of  prescribed  positive  scalar  weights  on  the  control,  u.  In  our  case,  R  and  u  are 
scalar  values. 

We  use  a  finite-volume  formulation  to  construct  a  finite-dimensional  approximation  to 
the  linear  dynamics  given  by  Eqn.  (2.32)  and  to  the  quadratic  spatial  integration  term  in  the 
cost-functional  (Eqn.  (2.39)).  The  variational  Hamiltonian  for  the  resulting  linear  quadratic 
control  problem  is  given  by 

H(  wN,  uN,  XN)  =  (/fQV/2  +  (nN)TRNnN/2  +  ( XNf  (A  V  +  B  V) ,  (2.40) 

where  XN  is  the  Lagrange  multiplier.  The  solution  to  the  steady  Riccati  equation  yields  a 
matrix  PN  which  maps  the  state  onto  the  Lagrange  multiplier  through 

XN  =  PNwN.  (2.41) 

The  steady  Riccati  equation  is  a  result  of  the  usual  necessary  optimality  conditions  for  this 
linear-quadratic  regulator  problem  (see  Paraskevopoulos  [12]).  For  notational  economy  at 
this  point  we  omit  the  superscript  N.  For  the  remainder  of  this  section  we  deal  only  with 
the  finite-dimensional  problem. 

The  Riccati  matrix  P  is  the  unique  positive-definite  solution  of 

PA  +  AtP  -  PBR_1BtP  +  Q  =  0.  (2.42) 

In  our  application,  the  Riccati  equation  is  solved  using  matlab,  and  the  control  is  given  by 

u  =  — R-1BrA  =  — R-1BtPw  =  -Kw.  (2.43) 

The  matrix,  K,  is  called  the  Kalman  matrix  or  the  gain  matrix. 


Feedback  control:  Numerical  results 

We  shall  consider  the  cost  functional  (Eqn.  (2.39))  with  the  following  weights  Qp(x)  = 
0,  Qu(x)  =  0.1,  Qp(x)  =  1  and  R  —  103.  These  weights  are  somewhat  arbitrary  but  reflect 
the  fact  that  we  want  to  limit  the  pressure  oscillations  and  that  in  the  open-loop  simulations 
the  velocity  responses  are  noticeably  smaller  than  the  pressure  responses.  In  this  low-Mach 
flow,  density  changes  are  expected  to  be  small. 

The  feedback  control  is  applied  to  linear  and  non-linear  systems.  The  pressure  response 
at  the  inflow  is  shown  in  Figure  2.44(a)  for  both  the  linear  and  non-linear  problems.  As  seen 
in  Figure  2.44(b),  the  control  required  to  damp  the  pressure  disturbances  is  much  smaller 
than  would  be  necessary  to  achieve  thermal  choking  (i.e.,  1000). 

The  same  control  is  applied  to  the  non-linear  system  at  time,  t  =  2.0  sec  after  the  limit 
cycle  has  been  established  (see  Figure  2.42(b)).  The  envelope  of  the  pressure  response  at 
the  inflow  is  shown  in  Figure  2.45.  The  oscillation  reaches  half-amplitude  in  approximately 
0.065  sec  with  an  initial  damping  rate  of  2.6  orders  per  second.  For  this  particular  finite- 
volume  representation,  the  control  is  very  effective;  however,  further  investigation  uncovers 
some  issues  of  concern. 
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Feedback  Control  Pressure  Response  to  a  1-D  Bombing  Pulse,  pout’=0 
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(a)  Inflow  pressure  response  to  feedback  control. 


Time  (sec) 


(b)  Feedback  control  history. 

Figure  2.44:  Linear  and  non-linear  feedback  control  response. 
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VACCG  Pressure  Response  with  Feedback  Control  at  t=2.0 

Ap/pj=  1  %,  y-1.396.  M-0.05, 2nd-Order  RK,  A t=10'5 


Figure  2.45:  Non-linear  response  to  feedback  control  applied  at  t  =  2.0 


Figure  2.46:  Pressure  gain  functional  at  different  mesh  densities. 
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Envelope  of  Pressure  Response  for  60-cell  Gains 
hp/pj=\%y  1.396,  Af=0.05, 2nd-0rder  RK,  A t=10'5 
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Figure  2.47:  Envelope  of  pressure  response  for  60-cell  gains  used  on  a  60-cell  mesh  and  a 
120-cell  mesh. 

Feedback  Control:  Reasons  for  concern 

Results  for  the  pressure  gain  functional  (KPN)  with  N  =  32,  60, 128  are  shown  in  Figure  2.46. 
Note  that  the  number  and  magnitude  of  the  oscillations  is  increasing  as  the  grid  is  refined. 
This  is  typical  of  weak  convergence  and  indicates  a  potential  difficulty  with  our  numerical  ap¬ 
proximation.  Specifically,  as  N  increases  the  finite-dimensional  operator  A  at  must  converge, 
in  an  appropriate  sense,  to  the  linear  approximation  of  the  continuum  PDE)  model. 
In  addition,  convergence  of  the  gains  in  the  LQR  problem,  also  requires  that  the  adjoints  of 
the  operators  converge.  In  other  application  of  distributed-parameter  control  theory,  lack  of 
adjoint  convergence  has  led  to  (only)  weak  convergence. of  the  functional  gains 

One  method  to  investigate  gain  convergence  is  to  use  the  60-cell  gains  on  a  finer  mesh 
of  120  cells.  That  is,  the  ith  gain  is  used  on  cell  numbers  2 i  and  2 i  +  1.  If  the  behavior  is 
the  same,  then  we  know  that  the  gain  is  grid-converged  -  the  converse  is  also  true.  A  plot 
of  the  pressure-response  envelope  is  given  in  Figure  2.47  which  shows  that  the  120-cell  case 
is  unstable  even  with  control. 

Distributed  Parameter  Systems 

In  the  previous  sections,  we  have  described  a  non-linear  finite-dimensional  (ODE)  model 
based  on  a  finite-volume  approximation  to  the  flow  equations.  The  non-linear  ODE  model 
was  linearized  and  used  to  construct  a  stabilizing  controller  based  on  an  linear-quadratic 
regulator  (LQR).  For  this  specific  discretization  the  state- feedback  controller  was  shown  to 
stabilize  the  non-linear  ODE  model.  However,  the  linear  controller  applied  to  a  non-linear 
ODE  model  based  on  a  finer  discretization’ did  not  result  in  a  stable  system.  Moreover, 
as  the  finite-dimensional  model  is  refined,  calculations  performed  on  the  sequence  of  linear 
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control  problems  indicate  that  certain  functional  gains  for  the  DPS  model  apparently  fail  to 
converge.  Thus,  we  are  interested  in  constructing  a  linear  distributed-parameter  model  and 
studying  control  and  approximation  issues. 


2.4.1  Units  and  Scaling 

Next,  we  introduce  non-dimensional  time  and  length  variables  given  by 

x 

t  =  y-j— ,  and  x  =  T, 

L/uq  L 


and  a  non-dimensional  mass  variable  given  by 


m  = 


m 


In  these  expressions  uo  and  po  are  convenient  values  of  flow  velocity  and  density.  The  scaled 
versions  of  the  (primitive)  dependent  variables  are 


p(t,  x )  =  p{tL/u o,  xL)/po 
u(t,x )  =  u(tL/uo,xL)/uo 
p(t,x)  =  p(tL/uo,xL)/(poul). 


The  linearized  PDE  for  the  scaled  primitive  variables  is 


d_ 

dt 


"  p ' 

u 

P  \ 

(t,  x)  + 

110 


[0  1/Mq  lj 


where 


d 

'P 

u 

(t,  x)  = 

dx 

. 

_P_ 

L 

0 

0 


L(T  —  l)Se(t,x)J 


se(t,  x)  = - $Se(iL/u0,xL). 

Pou0 


2.4.2  Boundary  Conditions-  -  ; 


(2.44) 


Classical  characteristic  analysis  of  Euler  flows  suggests  that  for  subsonic  flows  one  should 
provide  two  (Dirichlet)  conditions  at  the  inflow  and  one  at  the  outflow.  In  our  case  one  of 
the  inflow  conditions  is  used  to  model  unsteady  combustion  features  that  are  not  otherwise 
described. 


Inflow  Conditions 

At  x  =  0  the  mathematical  boundary  conditions  are: 

•  mass  flow  change  proportional  to  pressure  change 


d(pu)  dp  , 

o  fr  ®  o  ’ 

pUuU  pU 


•  constant  entropy  ( p(t$-t  =  K)- 
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In  terms  of  our  scaled  perturbation  variables,  the  first  condition  leads  to 

«((,  0)  +  Hi  0)  -  0)  =  0.  (2'45) 

while  the  second  condition  leads  to 

Mop(t,  0)  —  p(t,  0)  =  0.  (2.46) 

Outflow  Conditions 

At  the  outflow  x  =  1  the  only  boundary  condition  is 
•  fixed  pressure  ( p(t,L )  =  po), 
which  leads  to  the  simple  requirement 

p(i,  1)  =  0.  (2.47) 


2.5  Abstract  (DPS)  Model 


We  begin  with  the  Hilbert  space  Z  =  L2([0, 1],1R3),  which  consists  (essentially)  of  contin¬ 
uous  functions  on  the  unit  (spatial)  interval  with  values  in  1R3.  To  incorporate  boundary 
conditions  we  define  the  matrices 


Bl 

Br 


1  1  -a 

-1  0  Mq  J  ’ 

[0  0  1], 


(2.48) 

(2.49) 


where  “  =  (po/poWi)  (see  Eqn‘  2-45)‘ 

One  approach  is  to  build  the  boundary  conditions  into  the  domain  of  the  operator.  Define 
the  operator  A  :  D.  A )  Z  with  domain: 


D(A)  =  {tp  €  Z\ip  is  a.c.,ip'  6  L2,  BlV’( 0)  =  0  G  1R2,  BriP(  1)  =  0  €  IR} 


and  with  action  given  by: 


[M\  (*)  =  - 


1  1  0 
Oil 

Lo  1/M02  lj 


~l=zA 


']  (x). 


(2.50) 


The  homogeneous  case  of  Eqn.  (2.44)  can  then  be  written  as 

z(t)  =  Az(t), 


and  we  need  to  show  that  A  generates  a  Co  semigroup  on  Z. 

One  possible  approach  is  to  use  a  theorem  from  Pazy  [13,  see  Theorem  5.3,  page  20]. 
This  requires  that  we  show  that 
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1.  A  is  closed  and  densely  defined, 

2.  there  exists  real  numbers  M  >  0  and  u>  so  that  the  resolvent  set  includes  the  ray 
(u>,  oo),  and  the  resolvent  operator  R( A,  A)  satisfies  the  inequality 

for  $ft(A)  >  u ,  and  n  =  1,2,...  . 


2.5.1  Characterizing  the  spectrum  of  A 

Prom  Eqn.  (2.50)  we  have  [Atp](x)  =  -Aip'(x),  so  that  for  A  an  eigenvalue  of  A,  we  must 
have 

=  —\A~lip(x), 

where,  for  the  moment,  we  assume  that  A  is  invertible.  Prom  this  expression  we  have  the 
representation 

i/j(x)  =  exp  (— AA-1®)^^). 

The  boundary  conditions  (2.48,  2.49),  then  imply  that 

V»(0)  =  0  e  1R3.  (2.51) 

Hence,  eigenvalues  of  the  operator  A  are  equivalently  eigenvalues  of  the  matrix  M( A). 

To  construct  the  required  matrix  exponential  we  use  standard  results  from  one-dimensional 
Euler  flow  analysis  for  a  spectral  representation  of  the  matrix  A.  Specifically,  from  [14,  see 
pages  158,  ff\  the  eigenvalues  of  A  are 

<T  l  =  1,  <72  =  1  +  1/Mo,  cr3  =  1  —  1/Mo. 


Bh 

Hi?  exp  [— AA-1] 


Note,  in  particular,  that  for  0  <  M0  <  1,  the  matrix  A  is  well-defined  and  invertible. 
Additionally,  right  eigenvectors  («,)  of  A  are  given  by  the  columns  of  the  matrix  ^ 


1  M0/2  -Mo/2 

0  1/2  1/2 
0  1/(2M0)  -1/(2M0) 


while  the  reciprocal  left  eigenvectors  {yj)  are  given  by  the  rows  of 


'  1  0  -Mq  ' 
0  1  M0 
0  1  -M0 


Since  A  is  invertible,  it  follows  that  the  eigenvalues  of  A-1  are  the  reciprocals  of  the  cor¬ 
responding  eigenvalues  of  A  while  the  right  (left)  eigenvectors  of  the  two  matrices  are  the 
same.  Prom  these  observations  it  follows  that 

3 

exp  [-AA-1]  =  ^2  exP(— A/ (Tt)utv[ . 

t=l 
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Using  the  expressions  for  Br,  <7l,uu  and  vl ,  it  can  be  shown  that  the  last  row  of  the 
matrix  M( A)  in  Eqn  (2.51)  is 

[  0  M32(A)  M33(A)  ] , 

where, 

=  m.  (exp  (im) " cxp  (r^k))  - 

M33(A)  =  1  (exp  (f^r)  +  exp  (j^jjr)  )  • 

From  this  it  follows  that  the  characteristic  function  for  M (A)  is 

c(A)  =  (a  +  Mo-  Mq)  exp  -  (a- Mo  -  M02)  exp  •  (2‘52) 

2.5.2  Roots  of  the  Characteristic  Equation 
Define 

ci  =  a  +  Mo  —  Mq 

C2  =  ol  —  Mq  —  Mq  , 


and  note  that  ci  >  C2-  Additionally,  define 


Mo 

°  ~  1  +  Mo 


and 


_  Mo 
1-  M0' 


Recall  that  0  <  M0  <  1  so  that  b  >  a  >  0.  We  write  A  =  a  +  iu>  and  the  characteristic 
equation  (c(A)  =  0)  leads  to 


ci  exp(-ocr)  cos(aw)  -  C2  exp(fca)  cos(6w)  =  0  (2.53) 

ciexpf-ifcxlsinfa^)  ^.C2jxpfe).sih(i^)  =  ,P  •  j(2-54) 

Note  that  if  cos  (aw)  =  0  then  we  must  have  cos(6w)  =  0  so  that  b/a  is  rational.  We 
assume  that  this  does  not  occur  in  the  generic  case  ( i.e .  cos(aw)  cos(6w)  /  0).  From 
equations  (2.53,  2.54)  it  follows  that 


tan(aw)  +  tan(bw)  =  0. 


(2.55) 


Note  that  any  solution  of  the  characteristic  equation  is  a  solution  of  (2.55),  but  not  con¬ 
versely.  In  our  generic  case  we  use  the  trigonometric  identity 

tan(aw)  +  tan(fcw)  =  tan((a  +  6)w)(  1  -  tan(aw)  tan(6w)), 

where  we  have  assumed  that  cos((a  +  b)oS)  is  non-zero.  From  Eqn  (2.55)  it  follows  that 
tan(aw)  tan(6w)  <  0,  so  that  the  last  factor  is  is  not  zero,  which  implies  that  tan((a+6)w)  =  0, 
or  that 
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Eqn  (2.57) 

N  =  60 

N  =  128 

N  =  256 

CT/c 

Vk 

Ok 

Uk 

Ok 

.3905 

34.45 

.387 

34.1 

.392 

34.3 

.394 

34.4 

.3905 

103.35 

102.8 

.3905 

172.26 

171.4 

.389 

Bna 

.3905 

241.17 

.108 

238.5 

239.9 

.383 

gjggjl 

Table  2.3:  First  four  modes 


Next,  observe  that 


awfc  =  kn  —  boJk, 


so  that 

cos(awt)  =  (-l)*cos(M.)  => 

and,  similarly  that 

sin(acufc)  _  (_-,\k+i 
sin(6a;fc) 

Clearly,  cos  ((a  +  b)u>k)  is  non-zero.  With  u>k  values  from  (2.56)  we  formally  use  Eqns  (2.53, 
2.54)  to  solve  for  a.  The  results  are 


( 


=  log 


=  log  - 


C2  COS(6u/fc)  J 
clSm(au,k)\/(a  +  b) 


C2  sin(bcJk) , 


In  the  case  C1/C2  <  0  only  the  odd  values  of  k  yield  real  values  of  <7*,,  while  in  the  case 
if  C1/C2  >  0,  only  the  even  values  of  k  are  appropriate.  Finally,  if  ci  =  0  or  C2  =  0,  the 
spectrum  of  A  is.  empty.  With  Mq  >  0  the  case  =  &>,  —  0  cannot  occur.  In  summary,  we 


have 


log(|cx/c2)|  ±  i(2i  +  S) 7r 
(a  +  b) 


€  =  0,1,2,..., 


(2.57) 


where 


if  C1/C2  >  0 
if  C1/C2  <  0 


These  calculations  are  compared  to  values  from  the  the  LQR  model  for  the  case:  a  = 
.00288  and  Mo  =  .0455.  Note  that  b/a  «  2091/1909.  In  summary,  we  find 


\i  «  .124  ±  i34.45(2€  +  1),  €  =  0,1,... 

Table  (2.3)  presents  a  comparison  of  the  predictions  for  the  first  four  modes  from  Eqn  (2.57) 
and  from  the  LQR  methods  described  earlier. 
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2.5.3  A  Few  Interesting  Limits 

•  Neutral  stability  (9?(Ad  =  0)  occurs  when  |ci/c2|  =  1,  or  C2  =  —c\.  This  implies  that 
5  =  Mq  or  a  =  I/7. 

•  As  ci  |  0  we  have  -►  -00.  This  occurs  as  a  i  At  ci  =  0,  the  operator 

A  has  empty  spectrum. 

•  As  C2  t  0  we  have  5R(A^)  — >  00.  This  occurs  as  a  \  1^g  ■  At  C2  =  0,  the  operator  A 
has  empty  spectrum. 


2.5.4  Resolvent  of  A 

As  notes  at  the  beginning  of  the  Section  (see  page  59),  well-posedness  of  the  model  can 
be  established  through  certain  bounds  on  the  resolvent  operator.  Thus,  it  is  of  interest  to 
provide  a  representation  of  this  resolvent  operator.  We  have  performed  this  analysis  but 
omit  it  from  the  current  report. 


2.6  Boundary  Control 

The  physical  system  may  admit  the  possibility  of  control  through  the  flame  dynamics.  Our 
simplified  analysis  treats  the  fluid  as  a  perfect-gas,  while  the  flame  is  modelled  through  the 
inflow  boundary  conditions  (see  Section  2.4.2). 

Following  ideas  presented  in  [15]  we  generalize  the  mass  in-flow  boundary  condition  (2.45) 
to 

u(t,  0)  +  p(t,  0)  -  yp(t,  0)  -  yc(f,  0)  =  0,  (2.58) 

where  yp(t)  is  the  scalar  output  of  a  finite-dimensional  linear  system  driven  by  the  pressure 
perturbation  at  the  inflow  boundary,  and  yc(t )  is  the  scalar  output  of  a  finite-dimensional 
linear  system  driven  by  an  exogenous  control. 

For  the  pressure  driven  system  we  have 


fip(t)  =-Ar,up(t)±  ekp(t  ,Q)  _  _  ...  ...  ,  ......  ...  (2.59) 

yp(t)  =  Cppp{t)  +  (2.60) 

where  pp{i)  £  IRfc.  Without  loss  of  generality,  we  can  assume  the  system  is  in  control- 
canonical  form,  so  that,  for  example,  is  the  k  column  vector  with  zeros  except  for  unity 
in  the  last  entry.  The  transfer  function  for  the  system  (2.59,  2.60)  is 


Gp(s) 


YP(s) 

P(s) 


c^  +  cjs  +  ...  +  c£s(fc  x)  a 

dj,  +  djs  +  .. .  +  d£s(fc-1)  +  sk  (po/poUq) 


(2.61) 


Note  that  by  choosing  Cp  =  0,  j  =  1,2. . .  ,/c  we  recover  the  original  pressure  contribution 
from  Eqn  (2.45). 

Contributions  from  the  (new)  control  term  (v(t)  )  are  modelled  similarly 

pc(t)  =  Acpp{t)  +  etv{t) 
yc(t)  =  Ccpc(t), 


(2.62) 

(2.63) 
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where  /zc(f )  €  ffi/.  The  transfer  function  for  the  system  (2.62,  2.63)  is 

r  u\  -  =  Cc  +  Cc5  +  •  •  •  + 

c{)  V(a)  (11+4S  +  ...  +  disV-V  +  s* ' 

For  technical  reasons  Gc  is  strictly  proper;  there  is  no  direct  v{t)  »->  yc(t)  connection. 


(2.64) 


2.6.1  Second  Abstract  (DPS)  Model 

Define  the  Hilbert  space  W  =  IR^  x  ]Rfc  x  Z  with  elements  w  =  Up,  VO-  The  2nd  abstract 
model  is  given  by 

w(t)  -  Aww(t)  +  bwv(t),  (2.65) 

where  the  operator  Aw  :  D(AW )  t— >  W  has  domain: 

D(AW )  =  |/ic  €  IR e,fJ.p  €  JRk,xp  e  Z\ip  is  a.c .,xp'  €  £2, 

-  Ccnc  -  Cpy.v  +  ipi(0)  +  ^2(0)  -  5^3(0)  =  0, 

-^i(0)  +  M0V3(0)  =  0,  ^3(1)  =  0} 

and  action  given  by: 

fj-c  Acfic 

Aw  Up  =  Apfip  +  ekif>3(0)  ,  (2.66) 

xp  —Axp' 

where  the  3x3  matrix  A  has  been  defined  in  Eqn  (2.50).  The  control  operator  bw  :  IR  W 
is  given  by 

eev 

bwv  =  0  (2.67) 

0 

A  block  diagram  of  the  system  is  shown  in  Figure  2.48.  The  dark  vertical  bar  is  a  MUX 
operator;  the  output  vector  concatenates  the  input  vectors. 

2.6.2  Characterizing  the  spectrum  of  Aw 

From  the  definition  of  the  operator  Aw  in  Eqn  (2.66),  if  A  is  an  eigenvalue  we  must  have 

Acfj,c  =  A  fic  (2-68) 

ApUp  +  e*;V’3(0)  =  A  y,p  (2.69) 

-Aip'  =  Xxp.  (2.70) 

As  in  Section  2.5.1,  Eqn  (2.70)  implies  that 

xp(x)  =  exp(— XA~1x)xp(0). 

This  representation,  combined  with  the  boundary  conditions  and  Eqns  (2.68,  2.69)  leads  to 
-  XI -Ac  0  0  ]  /  He  \ 

0  A  I-Ap  e/t  [0  0 1]  (  ftp  I 


-Cp  =06  IRm+3, 

0  M(  A)  xp{  0) 

0  / 


(2.71) 


=M(  A) 
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To  Workspace 


Figure  2.48:  MATLAB  SlMULINK  Block  Diagram 

where  the  lower-right  3x3  block  (M( A))  has  been  defined  in  Eqn  (2.51).  We  are  interested 
in  characterizing  the  eigenvalues  of  the  matrix  (M( A))  on  the  left-side  of  Eqn  (2.71).  Note 
that  eigenvalues  of  Ac  that  are  distinct  from  the  eigenvalues  of  the  lower  right  (k+ 3)  x  (k-\- 3) 
block,  are  also  eigenvalues  of  M( A). 
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